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Abstract 

We introduce and analyze a phenomenological model for anti-angiogenic therapy in the treatment 
of metastatic cancers. It is a structured transport equation with a nonlocal boundary condition de- 
scribing the evolution of the density of metastasis that we analyze first at the continuous level. We 
present the numerical analysis of a lagrangian scheme based on the characteristics whose convergence 
establishes existence of solutions. Then we prove an error estimate and use the model to perform 
interesting simulations in view of clinical applications. 

Nous introduisons et analysons un modele phenomenologique pour les therapies anti-angiogeniques 
dans le traitement des cancers metastatiques. C'est une equation de transport structuree munie 
d'une condition aux limites non-locale qui decrit revolution de la densite de metastases. Au niveau 
continu, des estimations a priori prouvent Punicite. Nous presentons l'analyse numerique d'un schema 
lagrangien base sur les caracteristiques, dont la convergence nous permet d'etablir l'existence de 
solutions. Nous demontrons ensuite une estimation d'erreur et utilisons le modele pour produire des 
simulations interessantes au regard de possibles applications cliniques. 
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Introduction 



During the evolution of a cancer disease, a fundamental step for the tumor consists in provoking prolif- 
eration of the surrounding blood vessels and migration toward the tumour. This process, called tumoral 
neo-angiogenesis establishes a proper vascular network which ensures to the tumour supply of nutrients 
and allow the tumor to grow further than 2-3 mm diameter. It is also important in the metastatic process 
by making possible the spread of cancerous cells to the organism which then can develop in secondary 
tumors (metastases). Thus, an interesting therapeutic strategy first proposed by J. Folkman [15j in the 
seventies consists in blocking angiogenesis with the goal to starve the primary tumor by depriving it from 
nutrient supply. This can be achieved by inhibiting the action of the Vascular Endothelial Growth Factor 
molecule either with monoclonal antibodies or tyrosine kinase inhibitors. Although the concept of the 
therapy seems perfectly clear, the practical use of the anti-angiogenic (AA) drugs leaves various open 
questions regarding to the best temporal administration protocols. Indeed, AA treatments lead to rela- 
tively poor efficacy and can even provoke deleterious effects, especially on metastases [2JJ. Regarding to 
these therapeutic failures, it seems that the scheduling of the drug plays a major role. Indeed, as shown 
in the publication JT3], different schedules for the same drug can lead to completely different results. 
Moreover, AA drugs are never given in a monotherapy but always combined with cytotoxic agents (also 
named chemotherapy) which act directly on the cancerous cells. Again, the scheduling of the drugs seems 
to be highly relevant [23] and the optimal combination schedule between these two types of drugs is still 
a clinical open question. Thus, the complex dynamics of tumoral growth and metastatic evolution have 
to be taken into account in the design of temporal administration protocols for anti-cancerous drugs. 

In order to give answers to these questions, various mathematical models are being developed for 
tumoral growth including the angiogenic process. We can distinguish between two classes of models 
: mechanistic models (see for instance [5J [20]) try to integrate the whole biology of the processes and 
comprise a large number of parameters; on the other hand phenomcnological models aim to describe the 
tumoral growth without taking into account all the complexity levels (see [21] for a review and [TBI IT21 13] ) . 
Most of these models deal only with growth of the primary tumor but in 2000, Iwata et al. [TH] proposed 
a simple model for the evolution of the population of metastases, which was then further studied in [21 [5]. 
This model did not include the angiogenic process in the tumoral growth and thus could not integrate a 
description of the effect of an AA drug. We combined it with the tumoral model introduced by Hahnfeldt 
et al. 16J which takes into account for angiogenesis. The resulting partial differential equation is part of 
the so-called structured population dynamics (see [22] for an introduction to the theory) : it is a trans- 
port equation with a nonlocal boundary condition. Its mathematical analysis is not classical because the 
structuring variable is two-dimensional; as far as we know such models have only been studied in the case 
where one structuring variable is the age and thus has constant velocity (see [251 US])- This is not the 
case in our situation and the theoretical analysis of the model without treatment (autonomous case) was 
performed in [5]. 

In this paper, we present some mathematical and numerical analysis of the model in the non- 
autonomous case that is, integrating both cytotoxic and AA treatments and with a general growth 
field G satisfying the hypothesis that there exists a positive constant S such that G ■ v > S > where 
v is the normal to the boundary. We first simplify the problem by straightening the characteristics of 
the equation. We perform some theoretical analysis first at the continuous level (uniqueness and a pri- 
ori estimates) using the theory of renormalized solutions. Then we introduce an approximation scheme 
which follows the characteristics of the equation (lagrangian scheme). The introduction of such schemes 
in the area of size-structured population equations can be found in [JJ for one-dimensional models. Here, 
we go further in the lagrangian approach by doing the change of variables straightening the character- 
istics and discretizing the simple resulting equation, in the case of a general class of two-dimensional 
non-autonomous models. We prove existence of the weak solution to the continuous problem through the 
convergence of this scheme via discrete a priori L°° bounds and establish an error estimate in the case of 
more regular data. 

Finally, we use this scheme to perform various simulations demonstrating the possible utility of the 
model. First, as a predictive tool for the number of metastases in order to refine the existing classifica- 
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tions of cancers regarding to metastatic aggressiveness. Secondly, the model can be used to test various 
temporal administration protocols of A A drugs in monotherapy or combined with a cytotoxic agent. 



1 Model 

The model is based on the approach of [T51 [H E] to describe the evolution of a population of metastases 
represented by its density p(t, X) with X being the structuring variable, here two-dimensional X — (x, 9) 
with x the size (=number of cells) and 9 the so-called angiogenic capacity. It is a partial differential 
equation of transport type. The behavior of each individual of the population (metastasis), that is the 
growth rate G(t, X) of each tumor is taken from [TB] and is designed to take into account for the angiogenic 
process, as well as the effect of both anti-angiogenic (AA) and cytotoxic drugs (CT). Its expression will 
be established in the following subsection. The model writes 



d t p(t, X) + div(p(i, X)G(t, X)) = V(t, X) e]0, T[xQ 

(1.1) { -G-p(t,a)p(t,a)=N{a)f n l3{X)p(t,X)dX + f{t,a) V(t,tr) e]0,T[xdft 

p(0,X)=p°(X) VIefl. 

where f2, the birth rate P(X), the repartition along the boundary N(a) and the source term f(t,a) will 
be specified in the sequel, T is a positive time and v is the unit external normal vector to the boundary 
dft. 



1.1 The model of tumoral growth under angiogenic control (Hahnfeldt et al. 

US]) 

Let x(t) denote the size (number of cells) of a given tumor at time t. The growth of the tumor is modeled 
by a gompertzian growth rate modified by a death term describing the action of a CT. The equation is : 

dx f 9\ 

(1.2) — = gi(t,x) = aa; In ( - J - h'jc(t)H(x - x min ), 

where a is a parameter representing the velocity of the growth, 9 the carrying capacity of the environment, 
and the term h'-fc (t) H (x — x m i n ) stands for the effect of a cytotoxic drug, where 7c is the concentration of 
the CT, x m i n is a minimal size for the drug to be effective {x m i n > 1) arid the function H is a regularization 
of the Heaviside function (for example H(t) = l/2 + l/2tanh(i/i4T), with K being a parameter controlling 
the slope at 0), in order to avoid regularity issues in the analysis. The idea is now to take 9 as a variable 
of the time, representing the degree of vascularization of the tumor and called "angiogenic capacity". The 
variation rate for 9 derived in [TB] is : 

d9 

(1.3) — = g 2 (t,x,e) =cx - d9xi - er/ A (t)H(9 - 9 min ), 

dt 

where the terms cx and —d9x 2 ^ 3 represent respectively the endogenous stimulation and inhibition of the 
vasculature and ejA(t)H(9 — 9 m i n ) is the effect of an anti-angiogenic drug. The factor 2/3 comes from the 
analysis of [TCj which concluded that the ratio of the stimulation rate over the inhibition one should be 
homogeneous to the tumoral radius to the square. In the figure [TJ we present some numerical simulations 
of the phase plan of the system (|1.2[l . (|1.3[) . 

Following TB] , we assume a one compartmental pharmacokinetic for the AA and do the same for the CT 
(in [TB] there is no CT). We also assume that the drugs are administered as boli. This gives 

N 

7A W = ^z? Ae - c ^< t ^) J ff(t-0 

1=1 

where the tf are the administration times of the AA, Da is the administered dose and cWa the clearance. 
The expression for the CT is the same, with C instead of A. 
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Figure 1: Two phase plans of the system (jl.2p ~ p.3p . for different values of the parameters, without 

treatment (h = e = 0). In green, the nullclines. In both, b := (§) 2 = 17347. A. Parameters from [IB] : 
a = 0.192, c = 5.85, d = 8.73 x 1(T 3 . B. a = 0.192, c = 0.1, d = 1.4923 x 1(T 4 

1.2 Renewal equation for the density of metastasis 

We denote X = (x,6) and G{t,X) = ( gi (t, x, 9), g 2 (t, x, 9)). We define & = (f)^ and fi = (1,6) x (1,6) 
where b is the maximal reachable size and angiogenic capacity for (x,6) solving the system (jl.2p . (11.3p 
with initial size 1 (see [TT] for a study of this system without the CT term) . We consider that each tumor 
is a particle evolving in with the velocity G. Writing a balance law for the density p(t, X) we have 

dtp + div(pG) = 0, V(t,X) €]0,T[xQ 

that we endow with an initial condition p° e i°°(0). 

Metastasis do not only grow in size and angiogenic capacity, they are also able to emit new metastasis. 
We denote by b(cr, x, 9) the birth rate of new metastasis with size and angiogenic capacity a s dil by 
metastasis of size x and angiogenic capacity 9, and by fit, a) the term corresponding to metastasis 
produced by the primary tumor. Expressing the equality between the number of metastasis arriving in 
f2 per unit time (l.h.s in the following equality) and the total rate of new metastasis created by both the 
primary tumor and metastasis themselves (r.h.s.), we should have for all £ > 

(1.4) -/ p(t,a)G{t,a) ■ vda = [ [ b(a,X)p(t,X)dX + f(t,a)da. 

JdQ Jan Jn 

We assume that the emission rate of the primary and secondary tumors are equal and thus take /(t, a) = 
b(a,X p (t)) where X p (t) represents the primary tumor and solves the ODE system (|1.2D - (|1.3[) . We also 
assume that the new metastasis created have size x = 1 and that there is no metastasis of maximal size b 
nor maximal or minimal angiogenic capacity because they should come from metastasis outside of f2 since 
G points inward all along dft. An important feature of the model is to assume that the vasculature of the 
neo-metastasis is independent from the one which emitted it. This means that b(cr, X) — N(a)/3(x, 9) with 
N(<t) having its support in {a G dft; a = (1, 9), 1 < 9 < 6} and describing the angiogenic distribution 
of the metastasis at birth. We assume it to be uniformly centered around a mean value #0j thus we 
take N(l, 9) — 2Selee[e o -Ae,e o +A0], with A9 a parameter of dispersion of the new metastasis around #o- 
Following the modeling of [TH] for the colonization rate j3 we take 

P(x,6) = mi", 

with m the colonization coefficient and a the so-called fractal dimension of blood vessels infiltrating the 
tumor. The parameter a expresses the geometrical distribution of the vessels in the tumor. For example, 
if the vasculature is superficial then a is assigned to 2/3 thus making x a proportional to the area of the 
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surface of the tumor (assumed to be spheroidal). Else if the tumor is homogeneously vascularised, then 
a is supposed to be equal to 1. Assuming the equality of the integrands in (|1.4|) in order to have the 
equality of the integrals, we obtain the boundary condition of (jl.ljl . 



2 Analysis at the continuous level 

In the autonomous case, that is when G depends only on X and there is no treatment, the analysis of 
the equation has been performed in [5], It was proven the existence, uniqueness, regularity and 

asymptotic behavior of solutions. We present now some analysis on the equation (jl.ljl with a more general 
growth field G than the one defined in the section IT751 

Let n be a bounded domain in R 2 , with dft being piecewise C 1 except in a finite number of points. 
Let G : K x f2 — > M 2 be a C 1 vector field. We make the following assumption on G : 

(2.1) 3 5 > 0, G- u(t,a) > 8 > V < t < T, a e dfl. 
We do the following assumptions on the data : 

(2.2) p° g L°°(fi), P 6 L°°(fl), N e L°°(dSl), N > 0, f N{a)d<r = 1, / 6 L°°Q0, T[xdQ). 

Jan 

Remark 1. In the case of G being the one of the section \T7B if there is no treatment (that is, if e = h = 0, 
or t < t\) then we don't have G ■ u(t,a) > m > all along the boundary since G vanishes at the point 
(6,6). But since the problem was solved in this case (see f5jj) we consider that the time is the starting 
time of the treatment and that e or h is positive, which makes the assumption (j2.1|) true. 

Definition 1 (Weak solution). We say that p G L°°(]0, T[xf2) is a weak solution of the problem (jl.ip if 
for all test function <f> in C 1 ([0, T] x fl) with 4>(T, •) = 

/ p{t,X)[d t <f){t,X)+G{t,X)-'V<f>(t,X)]dXdt+ [ P °(X)(t)(0,X)dX 

(2.3) + / / {N(o-)B(t,p) + f(t,o-)}ct>(t,o-) = 

Jo Jon 

where we denoted B(t, p) := J n /3(X)p(t, X)dX . 

Remark 2. By approximating a Lipschitz function by C 1 ones, it is possible to prove that the definition 
of weak solutions would be equivalent with test functions in M^ 1,oo ([0,r] x Q) vanishing at time T. 

2.1 Change of variables 

Let X(t; r, a) be the solution of the differential equation 



(2.4) 



iX = G(t,X) 
X(t;t,ct) = a 



For each time t > 0, we define the entrance time r*(X) and entrance point cr t (X) for a point Xel! : 

t\X) := inf{0 < t < t; X(r;t,X) E Q}, a f (X) := X(r*(X); t, X). 
We consider the sets 

Cl\ = {!£!]; r*(X) > 0}, il\ = {X e f2; t\X) = 0} 
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lb 

Figure 2: The two changes of variables X\ and X 2 (represented only on the plane 0=1). 



and 



Q x := {(t,X) e [0,T] x!l;Ie Q 2 := {(t,X) e [0, T] x!l;Ie d 2 }. 



We also define Qi := {(*, r, cr); < r < t < T, a G 90} = X _1 (<3i) and notice that 

Ei := [0,T]x<9fi = {(t,X); t\X) = 0}, and S 2 = {(t, X(t; 0, a)); <t <T, a e dfl} = {(t,X); t\X) = 0}. 

See figure [5] for an illustration. We can now introduce the changes of variables that we will constantly 
use in the sequel. 



Proposition 1 (Change of variables). The maps 



(t,r,a) H- X(t;T,a) 



Q-2 



(t,Y) h-> X(t;0,Y) 

are bilipschitz. The inverse of X% is (t, X) 1— > (t, t (X), a (X)) and the inverse of X 2 is (t, X) M> (t, Y{X)) 
with Y(X) = X(0;t,X). Denoting Ji(t;r,a) = |det(DXi)| and J 2 {t;Y) = \dct(DX 2 )\, we have : 

/„ _n j / , \ i^, \ — \\ f* div G(u,X(u;T,a))du , T /, f* div G (u , X (u;0 ,Y)) du 

(2.5) Ji(t;T,<j) = \G(t,o-) ■ 17{o-)\e->T and J 2 (t; Y ) = e-io 

We refer to the appendix for the proof of this result and to the figure for an illustration. 
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Using these changes of variables we can write for a function / 6 i 1 (]0, T[xfi) 
f [ f(X)dX = f [ 

JO Jfl Jo Jo 



T r t 

f(X 1 (t;T,a))J 1 (t;T,a)dadr+ I I f{X 2 (t;0,Y))J 2 (t;Y)dY. 
dn JO Jfl 



We want to decompose the equation p.lj) into two subequations : one for the contribution of the boundary 
term and one for the contribution of the initial condition since they are "independent". Defining 

(2.6) Pi(i,r,a) :^p{t,X(t;T 7 a))J 1 (t;T,a) andp 2 (i,y) := p(t, X(t; 0, Y))J 2 (t; Y) 

we have, when the solution is regular : d{p\ — (ptP + div (pG)) J\ = and the same for p 2 . It is thus 
natural to introduce the following equations 



(2.7) 

where we denoted 



d t pi{t, T, a) = < t < t < T, a G dVt 

p 1 (T,T,*) = N(a)B(t,p 1 ,pa) + f(t 1 <T) < r < T, a e dQ 



and 
(2.8) 



B{t,p u p 2 ) = / / p(X(t;T,a))p l (t,T,a)dadT+ \ f3(X(t; 0, Y))p 2 (t, Y)dY, 
Jo Jan 



d t p 2 = t > 0, Y e n 

p 2 (0;Y) = p°(Y) Yen. 



We precise the definition of weak solutions to these equations. 

Definition 2. We say that a couple (pi,p 2 ) £ L°°{Qi) x L°°(]0, T[xCl) is a weak solution of the equations 
(|23 p -([2~% ]l if for all fa G C x (Qi) with fa(T, •) = we have : 

(2.9) III Mt,r,<j)dtfa{t,T,a)dadTdt+ [ [ \N(a)B(t,p 1 ,p 2 ) + f(t,a)}fa(t,t,a)=0, 
Jo Jo Jan Jo Jan 1 ' 

and for all fa £ C 1 ([0, T] x f2) with fa(T, •) = we have 

(2.10) f [ p 2 (t,Y)d t fa(t,Y)dt+ [ p Q (Y)fa(0,Y)dY = 0. 

Jo Jn Jn 

Remark 3. Ifpi is a regular function which solves (|2.7[) . then the weak formulation is satisfied since we 
have : 

/ / / pi(t,T,a)dtfa(t,r,a)dadTdt — / / fa(T,T,a)pi(T,T,a)dTda 
Jo Jo Jan Jo Jan 

V v ' 

=0 

T rt ~ _ f T f ~ 

t>x(t,T,a)dtpi(t,T,cr)dadTdt — / / fa(t,t,a)px(t,t,a)dadt. 

Jo Jan 

We prove now the following theorem, establishing the equivalence between the problem (II. lj) and the 
problem l|277 |) -l[2l5 |i . 

Theorem 1 (Equivalence between problem ((TTTJ) and problem (|^77 j) -(|2T5 ]> ). Let p 6 L°°(]0,T[xQ) &e 
a weafc solution of the equation (jl.ip . TTien (pi,p 2 ) given by (|2.6[) is a weafc solution of (|2.7[) - (I2.8[) . 
Conversely, if p\ and p 2 are weak solutions of (|2.7|> anc? (|2.8|l . t/ien £/ie function defined by 

(2.11) :=p 1 (i ) r*(X) )CT t (X))Jf 1 (t,r t (X), CT t (X))l Xens +p 2 (t,y(X))JJ 1 (t,y(X))l Xe ^ 
is a weafc solution of (jl.ip . 
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Proof. 

• Direct implication. Let p be a weak solution of the equation (jl.l)) . We will prove that /?2 defined 
by (EH) solves Let 4> 2 G C x ([0,T] x ft) with <fe(T,-) = 0. We define for left 2 (*,X) := 

<fi 2 (t,Y(X)) G W /1 ' 00 (Q2) and we intend to extend it in a Lipschitz function of [0,T] x ft so that we 
can use it as a test function in the weak formulation for p (see remark [5]). We define, for (i, r, a) G Qi, 
</>f(i,r, er) — (j) 2 (t, cr)C e (r) with C e (t) being a truncature function in C 1 ([0, +oo[) such that < £ e < 
1, £ £ (0) = 1, C e (r) = for r > e. Then 0f 6 and we set 0f(t,X) := 0f (i, t\X), cr*(X)) G 

W 1 ' 00 ((5i) since t*(X) and <x*(X) are Lipschitz from proposition [TJ We define then 

t>{ on Qi 
/>2 on Q 2 

The function </> e is Lipschitz on Q\, Lipschitz on Q 2 and <fr e G C([0, T] xft) since Qi<~}Q 2 — {{t, X); t 1 (X) = 
0}. Thus e G ^'^([O, T] x ft) with e (T, •) = 0. Using <j) £ as a test function in 1(23)1 . we have 



f p[d t <t>\ + G-V<t>\}dXdt+ [ 

jQi Jo 



nn 



{N(a)B(t, p) + f(t, a)} 0f (t, u)dtda 



+ [ p[dt<t>2 + G ■ V<p 2 ]dXdt + [ p°{X)(j) 2 (0,X)dX = = J* + I 2 . 
JQi Jn 

By doing the change of variables Xi in the term If and noticing that <f>\ (t, a) = 4>\(t, t, a) = <fi 2 (t, a)C e (t), 
we obtain 

I{= I I pi(t,T,a)d t 0i(t,a)( e (T)do-dTdt+ [ [ B(t, p)fc(t, a)? (t)da ► 0. 

Jo Jo Jo Jon 

Now doing the change of variables X 2 in the second term I 2 and noticing that dt4> 2 (t, Y) = dt(<fi 2 (t, X(t; 0, Y))) 
d t <fa(t, X(t; 0, Y)) + G(t, X(t; 0, Y)) ■ V0 2 (i, X(t; 0, Y)) gives the result. The equation on p t is proved in 
the same way. 

• Reverse implication. Let pi and p 2 be solutions of (|2.7p and (|2.8|) respectively. Define p(t, X) by 
dUnj, and consider a test function G ^([O, T] x ft) with 0(T, •) = 0. Then 0i := <p lQl G C^Qi), with 
0i (T, •) = 0, thus 4>i (t, t, er) := 0i(t, Xi (r, er)) is valid as a test function in the weak formulation of (|2.7p . 
In the same way 4> 2 (t,Y) := cj) 2 (t, X 2 (Y)) with cf> 2 := 0|q 2 is valid as a test function for (|2.8|) . Thus we 
have 



pi(t,T,a)dt4>i(t,T,a)dadTdt + / / B(t, pi, p 2 )(f>i{t,t 1 a)dadt 

Jo Jdn 



Qi 

+ I [ p 2 (t,y)dtMt,y)dtdy+ [ p°(y)fc(0,y)dy = Q 
Jo Jn Jn 

Doing the changes of variables gives the weak formulation of 1)1.1)1 . □ 

This theorem simplifies the structure of the problem (|1.1|) . In some sense, it formalizes the method of 
characteristics in the framework of weak solutions for our problem. The characteristics are straightened 
(see figure (2J) and the directional derivative along the field (t, G) is transformed in only a time derivative. 
Moreover, integrating the jacobians (which contains the transformation of areas) in the definitions of pi 
and p 2 , these functions are constant in time. The continuous analysis and discrete approximation of the 
problem is thus simplified. 

2.2 A priori continuous estimates and uniqueness 

In order to obtain a priori properties on the solutions of the equation, we will use the theory of renor- 
malized solutions first initiated by DiPerna-Lions [TU] in the case of K™ and further developed by Boyer 
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[7J in the case of a bounded domain. Let us first recall the result that we will use, which can be found in 
[7J. We need to introduce the following measure on ]0, T[xdSl : dpc = [G ■ v)dtda 

Proposition 2 (Renormalization property). Let p G L°°(]0, T[xf2) be a solution, in the distribution 
sense, to the equation : 

(2.12) dtp + div(pG) = 0. 

(i) The function p lies in C([0, T]; L p (fl)), for any 1 < p < oo. Furthermore, p is continuous in time 
with values in L°°(Cl) weak-*. 

(ii) There exists a function jp G L°°Q0,T[xdd;\dp,G\) such that for any h G C 1 (K), /or any test 
function tf> G C 1 ([0, T] x 51), and for any [to,ti] C [0,T], we /iaue 

' 1 / %)(fi t + div(Gty))dtdX + / h(p(t ))cl)(to)dX - f h{p{t x ))<t>{ti)dX 

/to JO -/O JO 

(2.13) - / / h{ 1P )(f>G ■ vdtda - [ [ h'(p)pdiv(G)<pdtdX = 

Jt Jao Jt Jo 

Remark 4. 

• By approximating the function s i— )• |s| 6j/ C 1 functions, it is possible to show that the formula (|2.13[) 
stands with h(s) — \s\. 

• TTie second point of the proposition implies in particular that h(p) has a trace which is h(~/p). 

• In this proposition is proved in the case of a much less regular field G but with the technical 
assumption that divG = 0, which is not the case here. Though, the proof can be extended to our case. 

Thanks to this result, we can prove the following proposition. 

Proposition 3 (Continuous a priori estimates). Let p G L°°(]0, T[xf2) be a weak solution of the equation 
(jl.ip . The following estimates stand 

(2.14) ||p(V)IUi(0) <eW"<-~\\ P \\ LHn) + fet-MW*-" f \f(s,a)\dads 

JO JdQ. 

and 

(2.15) IHlL-Q0,T[xn)<C oo 
with 

Coo = (|MU-||/9|U~|Mk- + II/IU-) ||G|k-e T ll« HvG l^- + \\p»\\ L ~e T W d ™G\\ L ~ 

Proof. 

• Estimate in L 1 . Let p be a weak solution of the equation (|1.1J) . Then in particular it solves (|2.12p 
in the sense of distributions. Thus the proposition O applies and gives a trace jp € L°°(]0, T[xdfl; {d^Q | ) . 
Now, by using (|2.13[) with h(s) = s and the definition of weak solutions to the equation (|l.lj) we have 
that for all G C}([0,T[xTi) 

[ -fp(t,a)4>(t,a)G(t,a) ■ vdodt = [ [ (jV(<r) f f3(X)p(t,X)dX + f(t,a)\ cf)(t,a)dadt 
io Jon Jo Jdn I Jo J 

which gives 

(2.16) - 1P {t,a)G{t,a) -v = N{a) [ (3{X)p{t, X)dX + f(t,a), a.e. 

Jo 
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In view of the remark 01 we know that \p\ is also a weak solution to the equation (|l.ip . with initial 
data \p°\ and boundary data \N(a)B(t, p) + f(t,a)\. By integrating this equation on O and using the 
divergence formula, we obtain in the distribution sense : 

4 f \p(t,X)\dX = - f G(t,a)-^p(t,a)\da= f \N(a)B(t, p) + f(t, a)\ da 
at Jn Jan. Jdn 

and thus 

j f \ P (t,X)\dX< U f \p(t,X)\dX + \f(t,a)\. 

A Gronwall lemma concludes. 

• Estimate in L°° . Using the proposition [2 we have pi and f>2 solving (12. 7\i and (|2.8p . By doing the 
changes of variables, using the definitions of pi and p2 and the formulas (|2.5p . we see that 

\\p(t, OlU^n) = IIp'iC*? •)IU 1 ao,t[xon) + l|P2(*> Ollz 1 ^)' v< > 

II II ^ ll~ II 11/^11 TlldivGlloo i |i~ II „T||divG||oo 

|p||i~(]0,T[xn) < ||Pl|| i0 o ( Q l) ll^l|L~(9n)e + ||P2||L-(]0,T[xO)e 

But solving explicitely the equation (|2.7p . we have 

|pi(*,T,o-)| = |?i(r,r,a)| - | N(a)B(t, pi, p 2 ) + /(t, a) 

< HJVlUlljSlloodlpiCr.Olk + llftCr.OIW + ll/IU- 

^lliviuii^iunp^oilix + ll/ll^ 



On the other hand, for p 2 we have ||p2|U°°Qo,T[xn) = 11/52(0) I U~(n) = \\p°\\L^(n)- □ 

Remark 5. TTie expression (|2.16p shows that in the case of a zero boundary data f , the trace 7/5 has 
some extra regularity, namely it is C([0, T}\ L l (d£l)). 

Corollary 1 (Uniqueness). If p and p' are two weak solutions of the problem (|1.1J) . then p = p' almost 
everywhere. 



3 Construction of approximated solutions and application to 
the existence 

In this section, we build a weak solution to the equation (|l.ip which, in view of the previous considerations, 
can be achieved by building a couple (pi,p 2 ) of solutions to the equations (|2.7p - (|2.8p (recall proposition 
[T]). We will achieve the existence by convergence of an approximation scheme to the problem (|2.7p - (|2.8p 
where the difficulty is restricted to the approximation of the boundary condition. Then we establish an 
error estimate in the case of more regular data. In order to avoid heavy notations, we forget about the 
tilda when referring to the problem (|2.7p ~ (|2.8p . We place ourselves in the case where = (l,b) 2 . 



3.1 Construction of approximated solutions of the problem (I2.7p - fl2.8p 

Let = to < ... < tk < ■■ < tK+i = T be a uniform subdivision of [0,T] with tk+i — tk = St. For 
the equation (|2.8p . we give ourself uniform subdivisions 1 = x\ < ... < xi < ... < Xl+i — b and 
1 = 6% < ... < 6 m < ... < 6l+i = b, with xi + i — xi = 9 m+1 — 9 m = 5x. The scheme for the equation (|2.8p 
is then given by : 

f P ° 2 (i, m) = ^ P °(x, e)d X de i<i, m <L 

\ pl +l {l,m) = pf(/,m) l<k<K,l<l,m<L 
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That is, p\(l 1 m) — p^ih m ) f° r au ^, h m - 

For the discretization of the equation (|2.7p , for each k let = tq < ... < Ti < ... < T k = t k with 

[0, 1] -> <9ft 

s 

,g € L 1 (cT2) we have J af ~ l g(a)da = f Q g(a(s))ds. Let = Si < ... < Sj < ... < Sm+i = 1 be an uniforr 
subdivision with s^+i — sj = 5a. The scheme is given by 



Ti+i — Ti = St. Let a : o a(s) ^ e & P arame ^ r ^ za ^ on °^ with |c'(s)| = 1 a.e., so that for 



(3.2) 



with 



and 
(3.3) 



p?((U) = N 3 B ((p° 2 ) ltm ) +/9 1 < j < M 

Pi +1 (h.i) = j) 1 < Jfe < A", < i < jfe, 1 < j < M 

P \ +1 {k + = N 3 B k +\p\ + \p k 2 +1 ) + 1 <j < M 



fc-1 M L 
i—1 j—1 l,rn — l 

~ f / /3(X(< fc ;r, ( r))p 1 (i fe ,r,a)dTda+ / /3(X(i fe ; 0, Y))p 2 {t k ,Y)dY 
Jo Jdn Jn 



Ph ■= sk IT IT P(X(t k ;r, a))dadr, ft m := (}(X(lk k ;0, (x, 6)))dxd8 

$ ■■= it i:r it f& ^ d ° dt > N i ■= h it N ^ d °- 



Notice that the schemes (13.11) and (|3.2J) are well-posed since the definition of Pi +1 (k + l,j) involves values 
of Pi +1 (i, j) only with < i < k. We denote by h — St + 8a + Sx and define now the piecewise constant 
functions px,h an d Pi.h on Q\ and [0, T[xf2 by, for < k < K , 1 < i < k, 1 < j < M and 1 < I, m < L 

p lyh (t,T,a(s)) = Pi(i,j) for t £ [t k ,t k+1 [, r efc-i.Ti], s £ [sj,s,-+i[ 
(3.4) pi ih (t,T,a(s)) = for t £ [t k ,t k+1 [, r e]t k ,t], s e [sj,s j+ x[ 

P2,h(t,X,6) = p\{l,m) for t £ [t k ,tk+i[,x £ [xi,xi+i[, £ [0 m ,0 m +i[. 

See the figure [3] for an illustration. Notice that we have 

t 




i=0 i=l 



Figure 3: Description of the discretization grid for Q\, only in the (r, t) plane. The arrows indicate the 
index used in assigning values to pi, ft in each mesh (formula (13.4p ). 
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k M M 

2 



(3.5) l|Pi,h(*fe:-)llii(]o,t fe [ x aa) = ^2^2\Pi( i d)\5t5a, ||/?2,h(*fc, Olk(«) = XI |p*(^»")| 

i—1 j—1 i,m=l 



Remark 6. 

• We take the same discretization step in x and 6 for pi but it would work the same with two different 
steps. 

• For more regular data, we could take point values instead of (|3.3J) . 

• It will be clear from the following that the scheme would converge the same regardless to the value 
that we give to p®(0,j). 



3.2 Discrete a priori estimates 

We prove the equivalent of the proposition [3] in the discrete case. Notice that there exists a con- 
stant C a such that T,jLi N 3 SfJ < I d n N ( a ) da + C " 5a = 1 + C ° 5a '= W N \\h and T,jLifj +1S(T < 
U°°ao,T[;i 1 (an)) + C a 5o- ■= \\f\\h- 



Proposition 4 (Discrete a priori estimates). Let (pi(i,j)) k and (p|(^ ?rl )) i . I being given by (|3.ip 
and (|3.2|) respectively. Then for all k 



\\P2,h{tk, OIU 1 ^) — Hp lk(«)> Hp2,ft||i°o(]o,T[xfi) - \\P IU°°(n) 
(3.6) ||p llh (t fc) Olk(]o,Mxau) < e^l^l-IWI" (||p°|k (n) + 



w\\L~\\N\\ h r 

(3.7) Wp^Wl^ < ll^lk-l|j8|U-™(llpi,h(**,-)lk i + llp IU i ) + ll/lli-- 

Moreover, if p° > then Pi(i,j), p2(l,Tn) > for all k,i,j,l,m. 

Proof. The non-negativity of the scheme is straightforward from the definition. The estimate for p2.h 
follows directly from the scheme (13. 2p . For the L 1 estimate on pi^ we compute, using the scheme (|3.1I) 

fe+1 M 

\\pi,h{tk+i,-)\\L^o,t h+1 [xdQ) =^2^2\p k+1 (i,j)\5tSa 

i=i j=i 

fe A/ M M 

= EEl^')l < M<7 + |^ fe+1 (^ +1 ,P2 +1 )| <s*$>;< jCT + **E fo 

1=1 j=l j=l j=l 

<||pi,/ l (tfe)||L 1 (]o, tfc [x 9 o) + | J B fc+1 ^ +1 ,P^ +1 )|^||A r IU + ^||/IU 

Now from the expression of B k+1 (p k+1 , p% +1 ) 

|s fc+i (^ +i ,pt +i )l< i i^i 1^11^(4,-) ik +\m\L-\\P2, h (t k ,-)\\n- 

Thus we obtain 

\\pi,h(t k +u 01k < (i + ll£lk°<Wk) ll/»i,/»(*fc, 01k + ll/3|k-<yt||^IUIIft ll h(*fc ) Oik + st\\f\\ h 

Now using a discrete Gronwall lemma we obtain 

II U Ml <- ll/3|U»l|JV|Lt fc /|| /. M| | H/3||L°°||A r |U||p2,/ I fa,Q||Li + \\f\\h \ 

\\pi,H(t k+ u-)\\^<e^ \\\pi*M\\» + \\P\\ L ~\\N\\ h / 
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Using ||pi,ft(to, -)Hli = and ||p 2 ,fe(*fc, Olli 1 ^) = II^IU 1 ^) ends tlie P r00r 01 tne 1/1 estimate. 
For the L°° estimate, we remark that 

IKhlli-ro,} -maxmax|pj(i,j)| = maxmax p^)iV, + / fe |) < ||7V|| L ^max |S fc (Pi, 

^ WD k i,j k j J k 

< ||^|U«||^|| £o oniax(||p lih (tfc,-)||Li + ||p2,fc(tfc,-)IU0 + 

k 



□ 



3.3 Application to existence of solutions to the continuous problem fl2.7l) - fl2.8|) 

Theorem 2 (Existence). Under the assumptions (|2.2p . there exists p\ € L°°(Qi) and p2 G L°°(]0, T[xil) 

smc/i £/ia£ pi /i — pi and p2 h P2 for the weak-* topology of L°° . Furthermore, (pi,p 2 ) is the unique 

h^O ' h-tO 

weak solution of (f277| ) -([2~%] ) ■ 

Proof. Uniqueness of the solution is straightforward for the problem (|2.8p and follows from the L 1 estimate 
on pi which can be derived following the proof of the proposition [3] The proof for the existence is rather 
classical and consists in passing to the limit in discrete weak formulations of (|2.7[) and (|2.8p . From the 
previous proposition, we obtain that the families {pi,/t} 5t s and {p2,h}g t § are bounded in L°° and 
thus there exist p\ e L°°(Qi), p2 G L°°(]Q, T[xSl) and some subsequences pi,/i„ and P2,h n such that 
Pi h„ - i Pi and p2 h n p2 for the weak-* topology of L°° . We have to prove now that (pi,p 2 ) 

is a weak solution of (|2.7p - (|2.8|) . The uniqueness of solutions to the equation implies then by standard 
argument that the whole sequence converges. It remains to prove that (/0i,p 2 ) solves (|2.7p ~ 
• The function p2 is a weak solution of (|2.8p . Let </> 2 be a test function for (|2.8p . We have 

ptfc+i r x i+i r^m+i 



f I P2, hn (t,Y)d t( f>2(t,Y)dYdt = J2J2 (4&m) f k+1 f 1+1 f m+1 dtMt^,e)dedxdt 

K L K L 



where we denoted $ 2 (£fc, Z, m) := ^jp- J^ +1 Jg m+1 2 (£fc, ^, 8)d0dx. Using the scheme (p 2 (Z, m) is constant 
in k) and < &2 (i_f^+i , I, rn) = since tn+i = T, we obtain 

/ / p2, hn (t,Y)dt<h{t,Y)dYdt = PzO 1 ' m )^2(T, I, m)(5x) 2 — £ p° 2 (l, m)* 2 (0, 1, m)(Sx) 2 

J ° Jn Lm=l l,m=l 

= -J2 P°2(l,m)<P2(0,l,m)(5xf = - [ p% hn (T)0(O, Y) — — > - / p°{Y)(/)(0, Y)dY 



l,m=l 

since /0 2 ft„ — ~ — * P°- Observing that the left hand side converges to Jq J fl P2dt<p2(t,Y)dYdt gives the 
result. 

• The function p\ is a weak solution of (|2.7p . Let <f>\ be a test function for (|2.7p . Then the same calcula- 
tion as above shows, with $>i(tk, ■= i f^ 3+1 <pi(tki T : cr)dadr and using that $>i(tx+i, = 
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as well as p k+1 (i,j) — p\{i,j) for 1 < i < k and 1 < j < M 

. K k M K k M 

L Pi,h n (t, t, a)d t (/)i (t, t, a)dadrdt = E E E Pl^' 3')*i(*k+l 5 h j)StSa - E E E PiihJl^k, i, j)Stda 

•'Q 1 k=l i=l j=l k=li=lj = l 

K M K-l k M 

= j2J2 P f(i,j)$i(t K+ i,i,j)6t5a+ Y,J2J2pi( i >ti® i ( tk + i > i >fi 5t5a 

i—l j—l k—l i—l j — l 

K-lk+1 M M 

k— 1 i—l j — l j — l 

K-l M M 

= -^^^+ 1 (fc + l,i)$l(« fe+ l,A; + l,i)^(T-^ / 9Kl.i)*l(*l 5 l:j)^ 
fc=l j=l j=l 
K M 

= - E E {N B\plpl) + /*) ^(t k , k,j)6tSa 
fe=U=i 

Defining the following piecewise constant functions : Bh(t, pi.h, Pi.h) = B k {p k , p\), Nh(a(s)) = Nj, fh(t,a(s)) - 
f k and ^i,h(t, ct(s)) = $>i(t k ,k,j) on [t k ,t k+ i[x[sj , Sj + i[, the previous equality reads 

pi thn (t,T,a)d t (f) 1 (t,T,a)dadTdt = / / (B hn (t, p 1Jln , p 2 ,h n )N hn (<r) + f hn (t,<j))$ 1 , hn (t,<j)dadt. 

Jst Jon 

We need the following lemma in order to conclude. 
Lemma 1. We have 

B hn {t,Pi,h n ,P2.h n ) ->■ B{t, Pl ,p 2 ) *-X°°(]0,T[). 

Proof. We define the piecewise constant function /3^(t, <t) as for JV), and and = /3f m for X G 

[. Let i G [t k ,t k+1 [, then 

M 

2 



B h (t, pi, h , P2,h) = B k (p k ,p k ) = / / (3l(T,cr)p hh (t,T,a)dTdcr-y^(3l J p k (k,j)dtScr+ V f3l m p k (l,m)(5 

J ° Jaa 3 =l i,m=l 

since we defined t, ct) = for r G]tfc, i]. Thus, for ip G -^ 1 (]0, T[) we have 

,T f T f* f K M f tk+1 

/ B h (t,p lth ,p 2 ,h)ip(t)dt= / / Pl{T,o)p 1 M{t,T,a)^{t)dadTdt-5tY,Y.PhPi( k ^) ^(t)dt5o 
Jo Jo Jo J an k=0 j =1 Jt k 

Pl(X)p2, h (t,X)il>(t)dXdt 



o Jo 

i} 

and we obtain the result by using p\ h n — 1 Pi *~L°° , p 2 h n — 1 P2*~L° c ,(3h n 5- | L8fc_ | < C 

and noticing that the second term goes to zero in view of the L°° bounds on p\ h (proposition [5]) and 
P. □ 

Using the lemma as well as Nu, fh n — 1 N, f * —L°°, \\Nh Jlx, 00 < C and $i h„ c ^ ' T ] xaa ^ M^t^Y 

hn—yo h n —>0 

the previous calculations give 

p lyhn (t,T,a)d t (f>i(t,T,a)dadTdt >- [ [ { N(a)B{t, p u p 2 ) + f(t,a)\ (f>(t,t,a)dadt. 

> l ^° Jo Jan 1 J 

On the other hand the left hand side also goes to f~ pi(t,r, a)dt<j>i(t,T, a)do-drdt. This proves that p\ 
__ 

verifies the definition [5] and ends the proof. □ 
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3.4 Error estimate 

We establish now an error estimate for the approximation of the equations (|2.7p - (|2.8p . For this section, 
we make the following assumptions on the data : 

(3.8) p° G W 1 ' 00 ^), /3 G W 1 ' 00 ^), N G W 1,o °(0fi), N>0, f N{a)dcr = 1, / G W^QO, T[x<9Q) 

It can be noticed that in order to perform the weak convergence of the approximated solutions and estab- 
lish theoretical existence to the continuous problem, we did not need to approximate the characteristics 
X(t;r,a) of the equation. In view of the error estimate though, we need to use another approxima- 
tion of the data than (|3.3p . For /3(X(t;r, cr)) we have to introduce an approximation Xh(t;r, a) of the 
characteristics given by a numerical integrator of the ODE system (|1.2l) - (|1.3p . Then we define 

'■= P(Xh(t k ; ruaj)), Pl m := (3(X h (t k ; 0, ( Xl ,9 m )) 
[ ' /* ■■= /(**,**), Nj:=N(aj). 

For gi and g% being two continuous functions on Qi and ]0,T[xf2 respectively, we define 

V 1 g 1 (t,r,a(s)) = gi(t k ,Ti,<Tj) for t G [t k ,t k+1 [, r G]ri_i,Tj], s G [s 3 , s j+1 [ 
Pigi(t,T,a(s)) = for t £ [t fc , r G]tfc,i], s G [s 3 , s 3+ i[ 

P 2 ff2(i, = g2(t k ,xi,6 m ) for t G [t fe ,£ fe+ i[, cc e]x ; ,x (+ i], 6 [6> m ,6» m+ i[ 

Lemma 2 (Projection error). Let (51,32) G M /1, ° (Qi) x H /1,oo (]0, T[xf2). TTiera there exists Cp 1 and 
Cp 2 such that 

(3.10) ||0i(tfc,-)-Piffi(tk,OIU-ao,t*D <C Pl /i, [|fl2(t*, — P 2 S2(**, 0IU-(n) <OpA 

We don't give the proof of this lemma since it is classical. We define e% t h := pi.h. — and e2.i1 '■= 
P2.f1 — P2P2 the errors of the schemes, with (pi,p2) solving the problem (|2.7p ~ (|2.8p . From the equation 
(lir7j) we have 

pi(tk+i,n,a-j) = pi(t k ,Ti,aj), < k < K, < i < k, 1 < j < M 
pi(ife + i,r fe+ i,crj) = N{uj)B(t k +i,px,P2) + f(tk,°j) 

and thus, subtracting this to (|3.1|) and denoting e k (i,j) — ei,fc(ifc, tj, £7j) we obtain 

/ ej +1 (i,j)-e{(i,j), 0<fc<X, 0<z<fc, 1<j<M 
[i ' U) { e k+1 {k + l, ] ) = N 3 E k +^+r k+1 

with 

fc M L 

£W = EE4-^ 1 ( i -i) Mff + E A 2 m 4 +i a,m)(fe) 2 

2—1 j — 1 /,m — 1 
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Nj {B k+1 ((pi(t k+1 ,Ti,a j )) iJ ,(p2(tk + i,x h 9 m )) l ^ - B(t k+1 , pi,p 2 )) 



approximation of the integral in B(t k , pi, P2) 



Hence the truncation error of the scheme comes only from the quadrature error coming from the 



Lemma 3 (Truncation error). Assume (j3T5) l and that (fioX^pi G ^'""(Qi), (/3oX 2 )p2 G W^QO, T[xO). 
TTien i/iere exists C r such that 

max |rf| < C r /i. 
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Proof. We have 

fc-l M L 

r ) = N AJ2 Z (#j " P{Xi(t k ;n, o-j))) Mh^ua^StSa + J2 (A 2 ™ - ^(t*; x h 6 m ))) 75 2 {tk,x h e m ) ( 

z— 1 j — 1 l,m= 1 

fc-l M L 

+ 5^ 5^ (tfc; 7-i, crj))J5i i r , a)StSa + f3(X 2 (t k ; xi,9 m ))p 2 (t kl x h 9 7n ) (5x) 2 

i—l j—1 l,m—l 

P(Xi(tk',T,cr))pt(tk,T,<T)dTdo- — I P(X 2 (tk,Y))p2{tk,Y)dY 

on 

f}(X\ (tk',r, a))pi(tk,T, a)drda}. 

/t fc _i Jon 

Thus 

Vj\ <||JV|| i »{|| i 8|| W r 1 ,oc -PiJTillL-HPxpilli! + ||X 2)h -P 2 X 2 || i oo||P 2 p 2 || il ) 

fc-l M 



+ Z / / |P2[(/3oX 2 )p 2 ](t fc ,r)-(/3oX 2 )p 2 (t fc ,x,0)|dx^+||(/3oX 1 )p 1 || LOO / 1 }. 

7 ~_ 1 J Xl J Xrr, 



Using the lemma [5] and the L 1 a priori estimate of proposition [3] gives the result. □ 

Remark 7 (Order of the truncation error). In order to have a better order for the truncation error we 
could use a more sophisticated quadrature method like for instance the trapezoid method on Q for y0 2 and 
on [0, tk-i[xd£l for pi (completed by a left rectangle method on [tk-i,tk[xd£l). Adapting the previous 
proof shows that if the numerical integrator used for the characteristics has order larger than 2, then the 
truncation error would have order 2 (order of the trapezoid method). 

Proposition 5 (Error estimate). Assume $£8% and that (pi,p 2 ) <E W 1,0O (Qi) x W^QO, T[xfi) is a 
regular solution of (|2.7[1 - (|2.8|) . Let p\^ and p 2 .h solve (13. 1|) and (|3.2p . Then there exists some constants 
C\ and C 2 such that 

( 3 - 12 ) \\pi,h{tk, •) - pi(tk, OIUMIOAD - C ± h > \\P2Ji(tk, •) - P2(tk, OlU^n) - 

Proof. In view of the lemma f2J it is sufficient to prove the proposition with W s p s (tk, ■) instead of p s (tk, •) 
(with s — 1,2). For the second estimate, we notice that 

\\P2M(t k , O-^frfa, -)\\ma) = ||e 2 ,fc(*fc, Oilmen) = E l e 2(^™)l (<^) 2 = E W> m ) ~ A^.M (5x) 2 

l.m i,m 

and the result follows from the definition of p 2 (/, to). For the first one, we have 

k M 

\\pi,h(tk-, ■) - PiPife, OIUmIo^D = IKh(*fci OlliHlo.UD = EE l e i(*>-?)| (5<(5cr - 

»=i j=i 

We can compute, using (|3.1ip 

ft M M M 

\\ei,h(tk+i,-)\W < EE \ e i +1 ^j)\ St5o- + \E k+1 \5t^Nj5o- + StY^ So- 

t=l 3=1 3=1 3=1 

< ||ei,h(tfer)l|i 1 + ^l|yS|U||^IU{||ei,h(tfe,0lli 1 + ||e2,/ l («fe+ir)|Ui} + aMt 

< (l + Jtll/SlloolliVllfcJIIei.fc^OllLi +C 2 U ||iV| \ h hP5t + C r h8t 

and conclude using a discrete Gronwall lemma. □ 
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Remark 8 (Order of the error). 

• By looking more carefully at the propagation of errors in the proof, we see that if we set p^l^m) = 
p°(l,m) (which is valid under (|3.8[) ). the error on p 2 comes only from the projection error. 

• If in addition, we follow the remark\T[for the approximation of the data, then the error between p\ t h 
and Pipi would be of order 2 if we had used a trapezoid method for the integral term in B(tk,pi,P2)- 

3.5 Application to the approximation of the problem (11. ip 

We explain now how we approximate the solution of fll.lj) from the approximation of the solutions of 
equations (|2.7p - (|2.8[) given by the schemes (|3.1[) . (|3.2[) . We translate formula (|2.11|) at the discrete level 
thanks to pij l: f>2,h given by (|3.4p and the solutions Pi(i,j), p^ihj) of the schemes (|3.ip and (|3.2p : 

(3.13) p h (t, X) p hh (t, t^X), a^X^J-Kt, r*(X), a*(X))l Xen t + p 2 ,h(t, K(X))^(t, Y(X))l Xen% . 

S ~* ' " ~v ' 

■ = Pl,h -=P2,h. 

The jacobians of the changes of variables Ji(t;r,a) — \G(t, a) -~^{a)\e^ dlvG ( u - x ( u > T > a ')) du an( j j 2 [t-,Y) = 

e I a dlvG ( u < x ( u '°> Y )) du are approximated respectively by J\ t h and J 2 ,hi piecewise constant functions con- 
structed similarly as in (|3.4p through Ji(i,j) :— e ri ^ k ^> and J% {l,m) := e 7~2(k,i,m) ^ wnere anc [ j~ 2 are 
one-dimensional quadrature methods such that 7i(k,i,j) ~ divG(X(s; r^, <jj))ds and T2(k,l,m) ~ 

J tfc divG(X(s; 0, (a;;, 9 m )))ds. The errors of these quadrature methods are denoted by r\, r^ and are 
assumed to be of order a%, 012 '■ 

n:=max \n(k,i,j)\ <C q (8t) a \ := max \r 2 (k,l,m)\ <C q (5t) a \ 

k,i,j ' k,l,m 

Hence, we have 

(3.14) J*(i,j) = Ji(* fcl r 4> ff i )e"'' l(M ' J ' ) , Jj(i.m) := e 7 " 2 ^^ = J 2 (f fe , z/, 0m)e-^ fe '<> m ). 

We define the following meshes : 

Vi(M,j) = {(t,X(t;r,a(s))); t G r e]r,-_i,Ti], s G [sj,s i+ i[} 

V 2 (k,l,m) = {(t,X(t;0,(xi,e m ))); t G [i*,tk+i[, as G G [9 m ,9 m+1 [} 

and, for a function g G C([0,T] x f2) 

(3.15) Fg(t,X) = g(t k ,X(t k ;Ti,o-j))l( tt x)eVi(k,i,j) + g(h, X{tk\ 0, (xi, 9 m )))l{t,x)ev 2 (k,i,m)- 

Remark 9. In the same way as the lemma there exists a constant Cp such that for all function 
g G W^°°{}0,T[xn) 

11.9 - Pfl'IU 1 (]o,T[xf2) < C P h. 

Theorem 3. Suppose that p G W '°°(]0, T[xQ) is a regular solution of the equation (jl.ip and Ze£ p/i 6e 
defined by (|3.13p . T/ien t/iere exists a constant C such that 

sup ||/9ft(t, •) -p(t,-)IUi(n) < C7i. 
te[o,T] 

Proof. In view of the remark^ it is again sufficient to prove the proposition with ¥p instead of p. Let t G 
[tfc,tfc+i[, then \\p h (t, -)-Pp(t, -)lki(n) = ||pi,/t(*fc, -)- F Pi(*k> Oll^n**) + ll/^.ft^fc: -Ppafafc, OIL^n**) 
with p s (t,X) := p(t, X)l X £ni (s = 1,2). We do the proof only for pi since it is similar for p 2 . We 
also don't write the dependency in a in order to avoid heavy notations. To obtain the complete proof it 
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suffices to add integrals with respect to a in the following and a in all the functions. Doing the change 
of variables X\ we have, noticing that Ppi(ifc, X(t k ; t)) = PiPi(ife, r)Pi Jf 1 (ifc, t) 



llPi.fcOfc,-) — Fpi(t A , -)ILi-(fj*fc) 



< 



|/0l,fc(*fc,T)| J ly l(tk,T)J 1 (t k ,T) - 1 

+ / |PiPi(t fe ,r)| ll-PiJf^tjb.rJJi^fc.rJldT 
Jo 

Now we have, using the definition (|3.14[) 



PlA t k,T)J-J l (t k ,T)-¥ 1 p 1 {t k ,T)P 1 J- 1 (t k ,T) J X {t k ,T)dl 



dr+ \pi,h(t k ,T) -¥ lPl (t k ,T)\dT+ 



— I IP T- 1 c - r i 



Jf e~^-/i - 1 < \e r 



\PJi 



Ui-PJil 



- 1 



Thus, since 



1| < 2rT, 



< oo from formula (|2.5|) and the fact that G ■ v > m > 0, and using |e 
there exists Cj such that 

ll^fc^fc.^^fc.T-)- 1 !!^ < C .A and||l-PiJ 1 - 1 (tj b ,T)Ji(i k ,r)||L- < Cj/i. 

The last inequality comes from the lemma [5] since Ji S W /1,1 ((0, +oo) ; L 1 (r)) 1, °° from the formula (|2.5p . 
Using then the continuous and discrete a priori L 1 estimates and the proposition \5\ gives the result. □ 

Remark 10. In the case of less regularity on the solution, we still have ph — 1 p, * — L°°(]0, T[xfi). 

Indeed, we write p h = pi.hJil+^^J^h = Jf 1 +P2,h J 2 X +Pi,h(Ji,h ~ J\)+P2,h{Ji,h ~ Jz)- Then we 



-l 
h 



as weH as 



< Ce 7 " 3 w'i/i C a constant. Using the theorem^ 



use that for s = 1,2 J, 
/or tte convergence ofpi^. and ~p2,h gives the result. 

Remark 11. In practical situations we are often only interested in the number of metastases and not in 
the density p itself. Notice that thanks to the formula J Q p(t, X)dX = J Q J gn P i(t, r, a)dadr+ J Q p°(X)dX , 
we don't have to compute the jacobians J%, J2 to get the number of metastases. Yet, we still have to 
compute the characteristics since they are requested in the computation of the boundary condition (see 
formula (|3.9p ). 



4 Numerical simulations 

4.1 Simulation technique and parameters 

Since our equation is two-dimensional, the computational cost can be relatively high because of the 
integral term in the boundary, especially for large-time simulations. In order to take into account that the 
metastases are born with a vasculature very close to a given value #0j we examine replacing the function 
N(a) by a dirac measure. In [BJ, we demonstrate that if we take N(l,9) = N £ (l,6) = ^leG[e -e .e +e] 
and let e go to zero the solution of the problem p.ip . in the case of an autonomous velocity field G 
and initial condition equal to zero, converges to the measure solution of a limit problem consisting in 
replacing N by a dirac measure in (1, 8q). We use here these considerations to reduce the computational 
cost and simulate only along the characteristic passing through (l,6>o), that is to say the scheme (|3.ip 
with only one discretization point o"o on dQ and N(o~) = l a=ao . Moreover, we use a Runge-Kutta method 
of order 4 for the approximation of the characteristics and a trapezoid method for the approximation of 
the boundary condition. 

The values of the parameters for the tumoral growth are taken from 16J, where they were fitted to mice 
data. Following [TH] and 0, we take a = 2/3 and fix the value of m arbitrarily. The values of the 
parameters (without the treatment) are gathered in the table [TJ The size (= volume) is expressed in 
mm 3 though until now it was thought as the number of cells. 
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a 

(day- 1 ) 


c 

(day- 1 ) 


d 

(day- 1 vol- 2 / 3 ) 


Xq (initial x) 
(vol) 


6 (initial 9) 
(vol) 


to 

(Nbof meta)(day- 1 )(vol- a ) 


a 


0.192 


5.85 


8.73 x 10- 3 


10~ B 


625 


10- 3 


2/3 



Tabic 1: Values of the parameters without treatment. 



4.2 Simulations without treatment 

A very important issue for clinicians is to determine the number of metastases which are not visible with 
medical imaging techniques (micro- metastases). Having a model for the density of metastases structured 
in size allows us to compute the number of visible and non-visible metastases. We took as threshold for 
a metastasis to be visible a size of 10 s cells, that is 100 mm 3 by using the conversion 1 mm 3 ~ 10 6 cells. 
In the figure [U we plotted the result of a simulation showing both the total number of metastases as well 
as only the visible ones. We observe that at day 20 the model predicts approximately one metastasis 
though it is not visible. At the end of the simulation, the total number of metastases is much bigger than 
the number of visible ones. 

Thus, an interesting application of the model would be to help designing a predictive tool for the total 
number of metastases present in the organism of the patient. In this perspective, we define a metastatic 
index as the integral of p on Q : 

MI(t) := [ p(t, X)dX. 
Jn 

Of course, this index depends on the values of the parameters, for example on the parameter to, as shown 
in the table [H The larger to, the larger the metastatic index. In this table, we remark that at least for 





MI(1.5) 


MI(7.5) 


MI(15) 


TO = 10~ 4 


5.80 x 10~ M 


6.60 x 10~ 2 


2.79 x lO" 1 


TO = 10~ 3 


5.80 x l(T a 


6.60 x lO" 1 


2.81 


TO = 10~ 2 


5.80 x 10- 1 


6.62 


30.1 



Table 2: Variation of the number of metastases with respect to to. 

times less than 15 days, it seems that the metastatic index is linear in to. Indeed, this can be explained 
by the fact that at the beginning, most of the metastases come from the primary tumor and not by the 
metastases themselves (see figure OA). This means that the renewal term in the boundary condition of 
(|1.1|) could be neglected for small times and that the solution of p.lj) is close to the one of 

dtp + div(pG) = 
-G-vp(t,a)=N(o)P(X p (t)) 
p°(X) = 0. 

But then, integrating the equation on Q gives MI(t) — J Q f3(X p (s))ds — to J q x p (s) a ds, where X p (s) = 
(x p (s), p (s)) represents the primary tumor and solves the system (|1.2p - (|1.3p with initial condition (xq, 9q). 
The figure 0B shows that for larger times metastases emitted by the metastases themselves are more 
important than the ones emitted by the primary tumor. The metastatic index for large time is then not 
anymore linear in to (result not shown). 

4.3 Simulations with treatment 

4.3.1 Anti-angiogenic drug alone 

We present various simulations of treatments, first involving an anti-angiogenic drug (AA) alone, 
in order to investigate the difference in effectiveness of various drugs regarding to their pharmacoki- 
netic/pharmacodynamic parameters. The first result shown in figure [S] takes the three drugs which were 
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5 1 1 5 20 25 30 35 40 



Figure 4: Evolution of the total number of metastases and of the number of visible metastases, that is 
whose size is bigger than 100 mm 3 (~ 10 s cells). 




B 



Figure 5: Number of metastases emitted by the primary tumor and by the metastases themselves. A. 
T=50. B. T=100 

used in [16] where only the effect on tumor growth was investigated, and simulates the effect on the metas- 
tases. The three drugs are TNP-470, endostatine and angiostatine and each drug is characterized by two 
parameters in the model : its efficacy e and its clearance rate cIta- These parameters were retrieved in 
[IB] by fitting the model to mice data. The administration protocols are the same for endostatine and 
angiostatine (20 mg every day) but for TNP-470 the drug is administered with a dose of 30 mg every 
two days. We observe that TNP-470 seems to have the poorest efficacy due to its large clearance. As 
noticed in [TB], the ratio e/clrA should govern the efficacy of the drug and its value is 0.13 for TNP-470 
and 0.39 for both endostatine and angiostatine. The model we developed is now able to simulate efficacy 
of the drugs on the metastatic evolution (figure [BJC). Interestingly, the drug which seems to be more 
efficient regarding to the tumor size at the end of the simulation (day 15), namely angiostatine, is not 
the one which gives the best result on the metastases. Indeed, the lower efficacy of endostatine regarding 
to ultimate size is due to a relatively high clearance provoking a quite fast rebound of the angiogenic 
capacity once the treatment stops. But since the tumor size was lower for longer time, the metastatic 
evolution was better contained. This shows that the model could be a helpful tool for the clinician since 
the response to a treatment can differ from the primary tumor to metastases, but the clinician has no 
data about micro- metastases which are not visible with imagery techniques. 

One of our main postulate in the treatment of cancer is that for a given drug, the effect can vary 
regarding to the temporal administration protocol of the drug, due to the combination of the pharma- 
cokinetic of the drug and the intrinsic dynamic of tumoral and metastatic growth. To investigate the 
effect of varying the administration schedule of the drug, we simulated various administration protocols 
for the same drug (endostatine). The results are presented in figure [7J We gave the same dose and the 
same number of administrations of the drug but either uniformly distributed during 10 days (endostatine 
2), concentrated in 5 days (endostatine 1) or in 2 days and a half (endostatine 3). We observe that the 
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tumor is better stabilized with a uniform administration of the drug (endostatine 2) but the number of 
metastases is better reduced with the intermediate protocol (endostatine 1). It is interesting to notice 
that again if we look at the effects at the end of the simulation, the results are different for the tumor size 
and for the metastases. The two protocols endostatine 1 and endostatine 2 give the same size at the end, 
but not the same number of metastases. Moreover, the best protocol regarding to minimization of the 
final number of metastases (endostatine 1) is neither the one which provoked the largest regression of the 
tumor during the treatment (endostatine 3) nor the one with the most stable tumor dynamic (endostatine 
2). In the figure [51 we investigate the influence of the AA dose (parameter Da) on tumoral, vascular and 
metastatic evolution. We observe that the model is consistent since it exhibits a monotonous response 
to variation of the dose. 

4.3.2 Combination of anti-angiogenic and cytotoxic drug 

An important question in clinical oncology is to determinate how to combine a cytotoxic drug (CT) that 
kills the proliferative cells and an anti-angiogenic (A A) drug which acts on the angiogenic process, either 
by blocking the angiogenic factors like VEGF (monoclonal antibodies, e.g. Bevacizumab) or by inhibiting 
the receptors to this molecule. The AA drugs are classified as part of the cytostatic drugs as they aim to 
stabilize the disease. For instance, in the treatment of breast cancer, patients which express the receptor 
HER receive a combination of Docetaxel (CT) and Herceptine (a tyrosine kinase inhibitor, AA). Two 
questions are still open : which drug should come before the other and then what is the best temporal 
repartition for each drug? Here, we perform a brief in silico study of the first question. Since we don't 
have real parameters for the cytotoxic drug we fix arbitrarily the value of each parameter h, clrc and 
Dc to 1, and perform simulations of the model to investigate combination of the CT and the AA. In the 
figure O we present the results of two simulations : one giving the AA before the CT (fig. OA) and the 
other one doing the opposite (fig. OB). Although in both cases the effect on the metastases is very good 
since the growth seems stopped (fig. E]D), it appears that the qualitative behaviors of the tumoral and 
metastatic responses are different regarding to the order of administration of the drugs (fig. [5]C and[5jD). 
According to the model, it would be better to administrate first the CT in order to reduce the tumor 
burden and then use the AA to stabilize the disease. Indeed, the number of metastasis at the end of the 
simulation is lower when the CT is applied first than in the opposite case. Of course, this conclusion 
depends on the tumoral growth and drugs parameters but this simulation shows that the model is able 
to exhibit different responses regarding to the order of administration between CT and AA drugs. 

5 Conclusion 

In this paper, we combined the models of fTB] and [TB] to obtain a model aiming at describing the effect 
of anti-angiogenic drugs on the metastatic growth. We established the well-posedness of the model and 
developed an efficient numerical scheme to perform simulations, which could be adapted to similar models 
in higher dimensions. The model can now be used in order to rationalize the temporal administration of 
the anti-angiogenic drugs. To achieve this, we have to implement the various pharmacokinetic models of 
the different AA drugs and then compare the in silico predictions to real patient data. 
An important open problem in this direction is the mathematical parameter identifiability of the model, 
that is to say the inverse problem of uniqueness of the parameters resulting in a given observation. It 
is also important to develop efficient numerical methods able to achieve the parameter identification 
from the data. Indeed, identifying the parameters m and a in a given patient could determine the 
metastatic aggressiveness of its cancer, through the metastatic index. This could lead to interesting 
clinical applications such as a refinement of the existing classifications like TNM or SBR, which deal only 
with the visible metastases. 

As shown in [T3], the metastatic response to AA treatment depends on the time schedule of the drug. 
The results of the simulations are encouraging in the perspective of using the model as a tool able 
to test various real temporal administration protocols of the drugs and to perform predictions of the 
mathematically optimized schedule for a given drug. Moreover, AA are never used in a monotherapy but 
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Figure 6: Effect of the three drugs from [IB]. The treatment is administered from days 5 to 10. Endo- 
statine (e = 0.66, cIta — 1.7) 20 mg every day, TNP-470 (e = 1.3, cIta — 10.1) 30 mg every two days 
and Angiostatine (e = 0.15, cWa = 0.38) 20 mg every day. A : Tumor size. B : Angiogenic capacity. C : 
Number of metastasis. 
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- Endostatine 1 


— e- 


- Endostatine 2 


— v- 


- Endostatine 3 



1.5 




Days 



c 



Figure 7: Three different temporal administration protocols for the same drug (Endostatine). Same 
dose (20 mg) and number of administrations (6) but more or less concentrated at the beginning of the 
treatment. Endostatine 1 : each day from day 5 to 10. Endostatine 2 : every two days from day 5 to 15. 
Endostatine 3 : twice a day from day 5 to 7.5. A : Tumor size. B : Angiogenic capacity. C : Number of 
metastasis. 
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Figure 9: Combination of an anti-angiogenic drug (AA) : endostatine (e = 0.66, cIta = 1.7 and Da — 
20mg) and a cytotoxic one (CT). The parameters for the CT are : h = 1, clrc = 1, and Dp = 
1. A. AA from day 5 to 10 then CT from day 10 to 15, every day. Tumor growth and angiogenic 
capacity. B. CT from day 5 to 10 then AA from day 10 to 15, every day. Tumor growth and angiogenic 
capacity. C : Comparison between both combinations on the tumor growth. D : Comparison between 
both combinations on the metastatic evolution. 



rather combined with a cytotoxic drug, and determining the best way to combine both drugs is still a 
clinical open question [23j . As shown in the figure^ the model could help in this direction, regarding both 
to tumor regression and metastatic evolution of the disease. We should also develop further the modeling 
in order to take into account for the competition effects between CT and AA. Indeed, by reducing the 
vasculature AA drugs should induce worse supply of both drugs and on the contrary some arguments 
are expressed in favor of a normalization effect on the tumor vasculature by AA therapy |19j , at least at 
the beginning of the treatment. These elements should be incorporated to the model via nonlinear terms 
involving the drugs concentrations in the equations Ql.2p - Ql.3p . The relative simplicity of the model (6 
parameters without treatment) is a great advantage in view of concrete applications since we have to be 
able to fit the model to patients' data in order to retrieve their parameters and then perform predictions 
about the optimal schedule. 

A fundamental problem that we have to integrate in our model is the one of toxicities which have to be 
dynamically controlled to optimize the scheduling of the drug. In the case of CT and on the tumoral 
growth, a model dealing with hematological toxicities is used to drive phase I clinical trials [2^1 [3] . In 
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our case, we also have to integrate a module to control the toxicity and address the resulting problem of 
optimization under constraints. 

Eventually, our model can be used to run in silico tests about the paradigm of metronomic chemotherapies 
which consists in delivering a cytotoxic drug at low doses and uniformly distributed in the treatment cycle 
rather than administrating the maximum tolerate dose (MTD) at the beginning of the cycle. Indeed, 
these metronomic protocols seem to have a dynamical anti-angiogenic effect [171 0] that can be integrated 
in the model of [16] for the tumour growth and in our model for the effect on metastases. 



A Proof of the proposition [T] 

The result for the second map is classical. For the first one, we have to deal with irregular points of the 
boundary 9f2. We denote by x the set of such points and set Xt '■— {X(t;T,^); £ £ X; < t < t}. In 
order to prove the result, it is sufficient to prove that the map 

t }o,t[xdci\ x -> n\\ Xt 

x " (t,<t) ^ X(t;T,a) 

is a diffeomorphism, that globally the map X[ : [0, t] x dfl — > £l\ is bilipschitz and that its inverse is 
X i-> (r'(X), ct'(X)). For the first point, since we avoid the irregular points of the boundary by excluding 
the set X; we have the C 1 regularity. It remains to prove that X[(t, a) is one-to-one and onto, and that 
its inverse is C . 

• The map X\ is one-to-one and onto. Let t > and X S H\. We have X\(t 1 (X) , a 1 (X)) — 
X(t; r*(X), a\X)) = X(t; T t (X),X(T t (X);t, X)) = X(t; t, X) = X. 

For the injectivity, we remark that if we have X(t;r,a) — X(t;r' .a') with for instance r' < r, then 
a = X(t; t', a') which is prohibited by the assumption that G ■ v{t, a) > 0. Thus X\ is one-to-one and we 
have, for (r, a) £ [0, t] x dfl : X(t; r*(Xf(r, a)), ct(X{(t, a))) = X(t; r, a) which implies t\X{{t, a)) = r. 
Thus, we have proven that the inverse of X\ is X H- (r*(X), a l (X)). 

• The map X\ is a diffeomorphism. We will prove the formula (|2.5J) for J\ which will conclude the 
proof by using the local inversion theorem. We have J\ (t; r, a) = \d T X\ A d a X\ | , with d a X[ := DyX o a' 
for a being a parametrization of dfl and DyX € Al2(K) the derivative in Y of X(t;r,Y) viewed as the 
flow on fL We compute 

d t {d T X{ A d a X{) = d T d t Xl A d a X{ + d T X\ A d t {.D Y X{ o a') = d T (G o X{) A d a X\ + d T X{ A DG o D Y X{ o ct' 
= DG o a r X^ A + d T X{ AflGo = div(G)(9 T ^ A d a X[). 

We compute now directly the value of Ji(t; t, a). We define 

T{K) = Xl^t + ^^-XKf^a) 

and now notice that we can write 



X{{t;t,a) = X\(t;t + h,X{{t + h\t,a)) 

= X((t;t + h,a) + D Y X{{t;t + h,o)(X\(t + h;t, a) - X{(t;t,o)) + o(h) 
= X[(t; t + h,cr) + hD Y X[(t; t + h,a)o G(t, a) + o(h). 

Now when h goes to zero D Y X{(t; t + h,a) — > D Y X{(t; t, a) — Id since X[(t; t, Y) = Y . Finally, we have 
T{h) -> -G(t,a), thus d T X[{t-t,a) = -G(t,a) and d T X\ A d a X\ (t; t, a) = -G(t,a)Aa' = G(t,a)-v(a). 
Solving the differential equation between times r and t and taking the absolute value then gives the 
formula l|Z5|) . 

• Globally, X\ is bilipschitz. It is possible to show that |||£>Wf |||i°°([o,t]x9f2) < e*'"' DG '" LOC (i ' T i xn . 
On the other hand, using the formula (DXl)^ 1 = Jf 1 t Com{DX[) and the fact that from (|2.5|) J^ 1 is 
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bounded on fl\ thanks to the assumption (|2.ip we have Mil roo < oo. Thus X\ and (XT) 1 

are Lipschitz on [0, t) x d£l\x and tt\ \xt respectively, and they are both globally continuous on [0, t] x 9f2 
and ft\. Hence they are globally Lipschitz. 

Remark 12. Using the same technique than in the previous proof, we can calculate the derivative of 
Xi(t; r, a) in the r direction. Indeed we compute, for all t, r, a 

Xi(t;r,a) = Xi(t;r + h,X x (T + h;r,a)) 

= Xi(t;r + h,a)+D Y Xx(t;T + h,a)(Xx(r + h;r, a) - X\(t; t, a)) + o(h) 
= X 1 (t; t + h, a) + hD Y Xi(t] r + h, a) o G(r, a) + o(h) 

which gives 

/ * -i \ avu \ t Xx(t;r + h,a) - Xx(t;r,a) 

(A.l) o T Xi(t; r, a) = hm = —DyXi(t; r, a) o G(r, a). 

h— >o h 
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